14. Análisis de Asociación de Variables Aleatorias

En esta sección, se presentan tests para analizar la asociación entre variables aleatorias continuas (ANOVA), o la asociación entre variables aleatorias discretas o categóricas (Chi-cuadrado, Bowker).

14.1. ¿Por qué ANOVA?

En los test de hipótesis estudiados en sesiones anteriores, consideramos el uso de la distribución \(t\)- student para analizar las medias muestrales. En lo que sigue utilizaremos el Análisis de la Varianza (ANOVA) o también llamado análisis de factores, para estudiar el efecto de uno o más factores (cada uno con dos o más niveles) sobre la media de una variable continua.

Test de medias para dos poblaciones normales con la misma varianza desconocida (caso común)

Sean \(X_1,\cdots,X_n\) e \(Y_1,\cdots,Y_m\) muestras independientes de poblaciones normales con medias desconocidas \(\mu_x\) y \(\mu_y\) y misma varianza desconocida \(\sigma^2\). Consideremos el test de hipótesis:

\[H_0: \mu_x = \mu_y\]

\[H_1: \mu_x \neq \mu_y\]

del Corolario del Teo de Fisher-Cochran se cumple:

\[ \frac{\bar{X} - \bar{Y} - (\mu_x - \mu_y)}{\sqrt{S_p^2(\frac{1}{n}+ \frac{1}{m})}} \sim t_{n+m-2}\]

donde

\[S_p^2 = \frac{(n-1)S_X^2 + (m-1)S_Y^2}{n+m-2}\]

de manera que se rechaza \(H_0\) si

\[\frac{\bar{X} - \bar{Y}}{\sqrt{S_p^2(\frac{1}{n}+ \frac{1}{m})}} > t_{\frac{\alpha}{2},n+m-2}\]

En el caso en que se requiera comparar más de 2 grupos, o examinar el efecto de 1, 2 o mas factores, este procedimiento se vuelve ineficiente, porque no queremos hacer un montón de t-tests para cada par.

Además, family-wise error rate (la tasa de error de la familia, error global) que es la probabilidad de cometer al menos un error de Tipo I en múltiples pruebas estadísticas realizadas en los mismos datos aumenta. Para \(c\) tests se calcula como:

\[\bar{\alpha} = 1-(1-\alpha)^c\]

• con c=2 tests, el error de tipo I es 0.0975 (alrededor de 2*0.05=0.1)

• con c=3 tests, el error de tipo I es 0.143 (alrededor de 3*0.05=0.15)

• con c=10 tests, el error de tipo I es 0.40 (alrededor de 10*0.05=0.5)

Por lo tanto, es mejor abordar con el modelo del Análisis de la Varianza (ANOVA).

14.2. Repaso

14.2.1. La distribución chi-cuadrado

Sean \(Z_1,\cdots, Z_n\, v.a. i.i.d. \, \sim {\it N}(0,1)\) entonces

\[ Y = Z_1^2+\cdots+Z_n^2 \sim \chi_{n}^2\]

donde \(n\) son los grados de libertad de la distribución.

Propiedades de la distribución \(\chi^2\):

(i) Propiedad aditiva: si \(X_1\) y \(X_2\) son dos v.a. independientes distribuidas \(\chi^2\) de \(n_1\) y \(n_2\) grados de libertad, entonces

\[X_1+X_2 \sim \chi_{n_1+n_2}^2\]

(ii) Esperanza y Varianza:

\[ E[X]= n, \qquad Var[X]= 2n\]

14.2.2. La distribución F

Sean \(X \sim \chi_n^2\) e \(Y \sim \chi_m^2\) dos v.a. independientes \(\chi^2\) de grados de libertad \(n\) y \(m\) respectivamente, entoncese se define:

\[F = \frac{\frac{X}{n}}{\frac{Y}{m}} \sim F_{n,m}\]

donde \(F_{n,m}\) es la distribución \(F\) de \(n\) y \(m\) grados de libertad. También se nota \(F(n,m)\).

suppressMessages(library(dplyr))
suppressMessages(library(plotly))
suppressMessages(library(ggplot2))
suppressMessages(library(rmarkdown))

vec <- seq(0,5,by=0.01)
params <- seq(1,20,by=1)
pvec <- list()
for (i in 1:length(params))
    for (j in 1:length(params)){
        k = length(params)*(i-1) + j
        pvec[[k]] <- df(vec,df1=params[i],df2=params[j],ncp=0)
}

steps1 <- list()
steps2 <- list()
fig <- plot_ly(width=600,height=600) %>% layout(title = "\n \n Densidad de Probabilidad F(n, m)",
                                                 yaxis = list(range=c(0,1)))
for (i in 1:length(params)){
     for (j in 1:length(params)){
        k = length(params)*(i-1) + j
        fig <- add_lines(fig, x=vec, y=pvec[[k]], 
                     visible=if ((i==1) && (j==1)) TRUE else FALSE,
                     mode='lines', line=list(color='blue'), showlegend=FALSE)
        steps2[[j]] = list(args = list('visible', rep(FALSE, length(params)*length(params))), 
                      label=params[j], method='restyle')
        steps2[[j]]$args[[2]][k] = TRUE
        steps1[[i]] = list(args = list('visible', rep(FALSE, length(params)*length(params))), 
                      label=params[i], method='restyle')
        steps1[[i]]$args[[2]][k] = TRUE
  }                   
}
fig <- fig %>% layout(sliders = 
                      list( list(active=0, currentvalue = list(prefix = "df1 (n): "), pad = list(t=20), steps=steps1),
                      list(active=0, currentvalue = list(prefix = "df2 (m): "), pad = list(t=100), steps=steps2)))
fig

Percentil

Sea \(F_{\alpha,n,m}\) tal que \(P\{F > F_{\alpha,n,m} \} = \alpha\), \(F_{\alpha,n,m}\) es el percentil \(100(1-\alpha)\) de la distribución \(F_{n,m}\).

Relación entre las distribuciones F y t

Si \(X \sim t_{n}\) entonces \(X^2 \sim F_{1,n}\). Por qué?

Por construcción, si \(X \sim t_{n}\), entonces

\[X = \frac{Z}{\sqrt{\frac{Y}{n}}}, \qquad con \, Z \sim {\cal N}(0,1) \,e \,Y \sim \chi_n^2\]

Entonces

\[X^2 = \frac{\frac{Z^2}{1}}{\frac{Y}{n}}, \qquad con \, Z^2 \sim \chi_1^2 \,e \,Y \sim \chi_n^2\]

14.2.3. La distribución de la media y varianza muestral del caso Normal

Teorema de Fisher-Cochran

Sean \(X_1,\cdots,X_n\) v.a. i.i.d. \({\cal N}(\mu,\sigma^2)\) entonces la media y varianza muestral cumplen:

\(\begin{equation} \begin{array}{lcll} (i) & \bar{X} &\sim& {\cal N}(\mu, \frac{\sigma^2}{n})\\ \\ (ii) & {\displaystyle \frac{(n-1)S^2}{\sigma^2}}& \sim& \chi_{n-1}^2 \\ \\ (iii)& \bar{X} &{\mathrel \perp} & S^2 \quad \text{(independentes)}\\ \end{array} \end{equation}\)

\(S^2\) es la varianza muestral:

\[S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i-\bar{X})^2 \]

Por qué en \( {\displaystyle \frac{(n-1)S^2}{\sigma^2}} \sim \chi_{n-1}^2\), el grado de liberdad es \(n-1\)?

14.3. ANOVA unidireccional (ANOVA de un factor)

Suponga que tenemos \(m\) tratamientos distintos, donde el resultado de aplicar el tratamiento \(i\) a un individuo es una v.a. normal. Estamos interesados en probar la hipótesis de que todos los tratamientos tienen el mismo efecto. Para ello se aplica cada tratamiento a una muestra distinta de \(n\) individuos y se analizan los resultados.

Formalmente, sean \(m\) muestras (grupos) independientes de tamaño \(n\):

\[X_{i1},\cdots,X_{in} \sim i.i.d. \,{\cal N}(\mu_i,\sigma^2), \qquad i=1,\cdots, m\]

Estamos interesados en probar:

\[H_0: \mu_1 = \mu_2 =\cdots =\mu_m\qquad\]

\[H_1: \exists i,j \in \{1,\cdots m \},\, \mu_i \neq \mu_j\]

14.3.1. Supuestos (y robustez)

Consideramos sus supuestos y robustez aquí:

• Normalidad: cada muestra (o v.a. \(X_{ij}\)) debe provenir de una distribución normal. Utilizar por ejemplo, test de Shapiro-Wilk. Utilizar Q-Q plot para descartar la presencia de outliers.

  • El ANOVA unidireccional se considera una prueba robusta contra el supuesto de normalidad, cuando el tamaño de cada muestra (o grupo) no es pequeño (y los tamaños de grupos son iguales); si no, consideremos test nonparamétrico como Kruskal-Wallis H Test.

• Varianza común (homogeneidad): misma varianza en los grupos. Se puede probar la igualdad de varianzas con el test de Levene o con el test de Brown-Forsythe.

  • Si este supuesto falla, consideraremos test nonparamétrico como Kruskal-Wallis H Test.

• Observaciones independientes. Este supuesto no puede ser probado con ningún estadístico, es una consideración de diseño

  • La violación de esto es muy grave (por ejemplo, la tasa de error de Tipo I está sustancialmente inflada).

Idea general de la derivación ANOVA

Se trata de construir dos estimadores de la varianza común, el primer estimador no depende de si \(H_0\) es cierto o no. En cambio, el segundo estimador asume que \(H_0\) es cierto, en caso contrario este estimador sobreestima a \(\sigma^2\). El test compara ambos estimadores y rechaza \(H_0\) cuando la tasa entre el segundo y primer estimador es suficientemente grande.

14.3.2. El primer estimador de \(\sigma^2\) que no depende de \(H_0\) (\(SS_W\))

Consideremos la suma de la cuadrados de las versiones estandarizadas de las v.a. involucradas:

\[\sum_{i=1}^m \sum_{j=1}^n \frac{(X_{ij} - \mu_i)^2}{\sigma^2} \sim \chi^2_{nm} \quad \text{ (Por qué?)}\]

Pero no conocemos \(\mu_i\), entonces podemos utilizar sus estimadores de máxima verosimilitud:

\[X_{i.} = \frac{\sum_{j=1}^n X_{ij}}{n}\]

Recordemos que:

\[ \sum_{j=1}^n \frac{(X_{ij} - X_{i.})^2}{\sigma^2} \sim \chi^2_{n-1} \quad \text{ (Por qué?)}\]

entonces por la propiedad de aditividad de la chi-cuadrado:

\[\sum_{i=1}^m \sum_{j=1}^n \frac{(X_{ij} - X_{i.})^2}{\sigma^2} \sim \chi^2_{nm-m}\]

Denominaremos

\[ SS_W = \sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{i.})^2\]

(\(SS_W\) es “within samples sum of squares”, o “residual sum of squares”)

Entonces

\[ \frac{E[SS_W]}{\sigma^2} = nm-m \quad \text{ (Por qué?)}\]

y luego

\[ \frac{E[SS_W]}{nm-m} = \sigma^2\]

con lo cual tenemos nuestro primer estimador insesgado de \(\sigma^2\):

\[MSE = \frac{SS_W}{nm-m} \quad \text{ (Por qué insesgado?)}\]

que se conoce como “error cuadrático medio” (mean squeare error), “cuadrado de la media de los residuos” (residual mean square) o “varianza intragrupos” (within-group variance). Notar que no considera ningun supuesto sobre \(H_0\).

14.3.3. El segundo estimador de \(\sigma^2\) que depende de \(H_0\) (\(SS_B\))

En este caso suponemos que se cumple \(H_0\), es decir

\[\mu = \mu_1 = \mu_2 = \cdots = \mu_m\]

De manera de que las medias muestrales cumplen:

\[X_{1.},X_{2.},\cdots,X_{m.} \sim i.i.d. \, {\cal N}\left(\mu, \frac{\sigma^2}{n} \right) \quad \text{ (Por qué?)}\]

de donde se pueden derivar distribuciones normales estandarizadas

\[Z_i = \frac{X_{i.}-\mu}{\sqrt{\frac{\sigma^2}{n}}} \sim {\cal N}(0,1) \qquad i=1,\cdots,m\]

de manera que:

\[\sum_{i=1}^m Z_i^2 = \sum_{i=1}^m \frac{(X_{i.}-\mu)^2}{\frac{\sigma^2}{n}} \sim \chi^2_{m}\]

Pero no conocemos \(\mu\). Si lo reemplazamos por su estimador de máxima verosimilitud:

\[X_{..} = \frac{\sum_{i=1}^m X_{i.}}{m} = \frac{\sum_{i=1}^m \sum_{j=1}^n X_{ij}}{nm}\]

entonces

\[\sum_{i=1}^m Z_i^2 = \sum_{i=1}^m \frac{(X_{i.}-X_{..})^2}{\frac{\sigma^2}{n}} \sim \chi^2_{m-1} \quad \text{ (Por qué?)}\]

Denominaremos

\[ SS_B = n\sum_{i=1}^m (X_{i.} - X_{..})^2 = \sum_{i=1}^m \sum_{i=1}^n (X_{i.} - X_{..})^2\]

(\(SS_B\) es “between samples sum of squares”, o “model sum of squares”.)

Entonces, si \(H_0\) es cierto, se cumple

\[ \frac{E[SS_B]}{\sigma^2} = m-1 \quad \text{ (Por qué?)}\]

y luego

\[ \frac{E[SS_B]}{m-1} = \sigma^2\]

con lo cual tenemos nuestro segundo estimador insesgado de \(\sigma^2\):

\[MSB = \frac{SS_B}{m-1} \quad \text{ (Por qué insesgado?)}\]

que se conoce como “cuadrado de la media entre muestras” (between samples mean square), “cuadrado de la media del modelo” (model mean square) o “varianza entregrupos” (between-group variance).

14.3.4. Test estadístico

Se define como estadístico:

\[TS = \frac{MSB}{MSE} = \frac{varianza \, entregrupos}{varianza \, intragrupos}\]

Cuándo rechazamos la hipótesis nula?

Se rechazará la hipótesis nula cuando TS sea suficientemente grande.

¿Cuál es la distribución de TS?

\[ TS \sim F_{m-1,nm-m}\]

Demostración:

\[ TS = \frac{MSB}{MSE} = \frac{\frac{SS_B}{(m-1)}}{\frac{SS_W}{(nm-m)}} = \frac{\frac{n}{m-1}\sum_{i=1}^m (X_{i.} - X_{..})^2}{\frac{1}{nm-m}\sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{i.})^2} = \frac{\frac{1}{(m-1)}\sum_{i=1}^m \frac{(X_{i.} - X_{..})^2}{\frac{\sigma^2}{n}}}{\frac{1}{(nm-m)}\sum_{i=1}^m \sum_{j=1}^n \frac{ (X_{ij} - X_{i.})^2}{\sigma^2}} \]

Sea entonces \(F\sim F_{m-1,nm-m}\) y \(F_{m-1,nm-m,\alpha}\) tal que:

\[P\{F > F_{m-1,nm-m,\alpha}\} = \alpha\]

Se calcula entonces el valor de \(TS\), si \(TS > F_{m-1,nm-m,\alpha}\) se rechaza \(H_0\).

14.3.5. Identidad de suma de cuadrados

Definimos \(SS_{total} = SS_{T}\) (total sum of squares) como:

\[\sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{..})^2\]

Cumple:

\[\sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{..})^2 = \sum_{i=1}^m \sum_{j=1}^n (X_{ij} - X_{i.})^2 + \sum_{i=1}^m \sum_{j=1}^n (X_{i.} - X_{..})^2 \quad \text{ (Por qué?)}\]

Así que tenemos:

\[SS_{total} = SS_{intragrupos} + SS_{entregrupos}\]

\[SS_{T} = SS_{W} + SS_{B}\]

y además

\[\sum_{i=1}^m \sum_{j=1}^n X_{ij}^2 = nmX_{..}^2 + SS_B + SS_W \quad \text{ (Por qué?)}\]

Ejemplo (industria de automóviles)

Una industria de automóviles desea probar la calidad de 3 tipos de combustibles. Para ello observa el rendimiento (millas por 10 galones de combustible) de los 3 tipos de combustibles en 15 vehículos idénticos (5 con cada tipo de combustible). (Asumimos que los supuestos de ANOVA se mantienen aquí, pero hay que examinar los supuestos en tu análisis de los datos.)

Cuál es la hipótesis nula aquí?

Veamos primero la distribución de los datos:

library(tidyr)
library(ggplot2)
library(dplyr)
library(cowplot)

options(repr.plot.width=10, repr.plot.height=5)
options(digits=4)

gas1 <- c(220,251,226,246,260)
gas2 <- c(244,235,232,242,225)
gas3 <- c(252,272,250,238,256)

df <- do.call(rbind, Map(data.frame, gas1=gas1, gas2=gas2, gas3=gas3))
df <- df %>% pivot_longer(cols=c('gas1', 'gas2', 'gas3'),names_to='combustible', values_to='milla')
df %>% group_by(combustible) %>% summarise(M=mean(milla), SD=sd(milla))

p1 <- ggplot(df, aes(x=combustible, y=milla)) + geom_violin(trim=FALSE) + 
    geom_dotplot(binaxis='y', stackdir='center', dotsize=1, binwidth=1.5) +
    stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1), geom="pointrange", color="red", lwd=1) + 
    theme(text=element_text(size=14))
p2 <- ggplot(df, aes(x=combustible, y=milla)) + geom_boxplot() + theme(text=element_text(size=14))

plot_grid(p1, p2, align="h", nrow=1, ncol=2, rel_widths= c(1/2, 1/2))

A tibble: 3 x 3

combustible	M	SD
<chr>	<dbl>	<dbl>
gas1	240.6	16.965
gas2	235.6	7.701
gas3	253.6	12.280

En tu tarea, tienes que desarrollar el test. Aquí, mostramos el resultado de R:

fm <- aov(milla ~ combustible, data=df)
summary(fm)

            Df Sum Sq Mean Sq F value Pr(>F)
combustible  2    863     432     2.6   0.12
Residuals   12   1992     166               

¿Cómo calculamos el p-value arriba con la función pf en R?

¿Cómo calculamos el valor crítico (\(\alpha=0.05\)) con la función qf en R?

14.3.6. Caso de muestras con diferentes tamaños

¿Qué ocurre si las muestras no son del mismo tamaño?

\[\sum_{i=1}^m \sum_{j=1}^{n_i} \frac{(X_{ij} - X_{i.})^2}{\sigma^2} \sim \chi^2_{N} \quad \text{ (Por qué?)}\]

donde \(N = \sum_{i=1}^m n_i - m\). Y entonces,

\[ SS_W = \sum_{i=1}^m \sum_{j=1}^{n_i} (X_{ij} - X_{i.})^2\]

por otra parte

\[\sum_{i=1}^m \frac{(X_{i.}-X_{..})^2}{\frac{\sigma^2}{n_i}} \sim \chi^2_{m-1}\]

y se definen

\[ SS_B = \sum_{i=1}^m n_i(X_{i.}-X_{..})^2 = \sum_{i=1}^m \sum_{j=1}^{n_i} (X_{i.}-X_{..})^2\]

\[ TS = \frac{\frac{SS_B}{m-1}}{\frac{SS_W}{N}} \sim F_{m-1,N} \qquad (N= \sum_{i=1}^m n_i - m)\]

Finalmente, si \(TS > F_{\alpha,m-1,N}\) se rechaza \(H_0\).

14.4. Métodos de comparaciones múltiples

Si el test ANOVA es estadísticamente significativo, lo único que podemos concluir es que una o mas de las diferencias entre grupos es significativa, pero no sabemos que grupos son los que difirien. Para ello es necesario realizar “Comparaciones múltiples”. Se trata de comparar cada grupo con otro grupo, o un promedio de grupos se puede comparar con otros. Se consideran dos posibilidades:

• Comparaciones planeadas: existe interés por algunos grupos en particular

• Comparaciones post hoc: no existen hipótesis específicas de algunos grupos

Algunas preguntas

Si el test ANOVA es estadísticamente significativo, implica que ¿existe al menos un grupo cuyo comportamiento difiere de los otros?

No necesariamente, podría ocurrir que sea la combinación de grupos que provoque las diferencias. El test de Scheffe permite pesquisar las diferencias subyacentes.

¿Es posible encontrar diferencias siginificativas con las comparaciones múltiples, aunque ANOVA no haya sido significativo?

Si, es posible, debido a que los tests de comparaciones múltiples están mas enfocados (tienen más potencia) que ANOVA.

¿Es útil el resultado del test ANOVA?

Claro, cuando corresponde a la hipótesis en estudio!

A menudo, sus preguntas experimentales están más enfocadas en algunos grupos. En estos casos, puede saltar directamente a los tests de comparaciones múltiples (aunque algunas personas no están de acuerdo con esta afirmación). Tenga en cuenta que todos los cálculos de comparación múltiple utilizan MSE de la tabla ANOVA. Por lo tanto, incluso si no le importa el valor de \(F\) o el p-value en ANOVA, los tests de comparación múltiple aún requieren que se calcule MSE en la tabla ANOVA:

\[MSE = \frac{SS_W}{\sum_{i=1}^m n_i - m} = \frac{\sum_{i=1}^m \sum_{j=1}^{n_i} (X_{ij} - X_{i.})^2}{\sum_{i=1}^m n_i - m}\]

Aquí presentamos algunos métodos de comparaciones múltiples. Tukey y Bonferroni son los más comúnmente usados, pero LSD es más fácil y ayuda a entenderlos.

14.4.1. LSD

En diferencias menos significativas, o least significant Differences (LSD), se define el siguiente estadístico cuya distribución es t-student, bajo

\[H_0: \mu_1 = \mu_2 = \cdots = \mu_m\]

\[H_1: \exists i,j, \mu_i \neq \mu_j\]

Para cada par de grupos i,j se calcula

\[t_N = \frac{X_{i.} - X_{j.}}{\sqrt{MSE \left(\frac{1}{n_i}+ \frac{1}{n_j}\right)}}, \qquad df = N = \sum_{i=1}^m n_i - m\]

Si \(|t_N| > t_{\frac{\alpha}{2},N}\) se rechaza \(H_0\) en favor de \(H_1: \mu_i \neq \mu_j\).

O calculamos el p-value:

\[p_{LSD} = 2P(T_N > |t_N|)\]

Si \(p_{LSD} < \alpha\) se rechaza \(H_0\) en favor de \(H_1: \mu_i \neq \mu_j\).

Aquí \(\alpha\) es igual que \(\alpha\) original (e.g., 0.05), y no se realiza ningún ajuste para el número de comparaciones. Este enfoque presenta problemas con el error de tipo I, en función del número de comparaciones posibles.

14.4.2. Bonferroni

Utiliza el mismo estadístico que LSD, pero calcula la significancia \(\alpha\) considerando que se realizan c comparaciones.

Dos posibilidades:

• Se corrige el p-value multiplicando por el número de comparaciones, y lo compara con el \(\alpha\) original (e.g., 0.05):

\[p_{Bonferroni} = c* p_{LSD} = \frac{m(m-1)}{2} *p_{LSD}\]

• O se calcula un nuevo valor de \(\alpha\) para mantener el error de tipo I global:

\[ \alpha_{nuevo} = \frac{\alpha}{c} = \frac{\alpha}{\frac{m(m-1)}{2}}\]

the Por ejemplo, si \(m=4, c=6\) y \(\alpha=0.05\), entonces \(\alpha_{nuevo} = 0.0083\)

Bonferroni puede ser demasiado conservador (\(p_{Bonferroni}\) bastante grande, o \(\alpha_{nuevo}\) bastante pequeño), lo que dificulta detectar diferencias que realmente existen (reduciendo su potencia). Pero puede garantizar el error de Tipo I, y es facíl entender la corrección. También lo usamos en comparaciones planeadas.

14.4.3. Tukey y Tukey-Kramer

Tukey: mismo tamaño muestral n

Se define el siguiente estadístico:

\[q_N =\frac{X_{i.} - X_{j.}}{\sqrt{\frac{MSE}{N}}}, \qquad df = N = \sum_{i=1}^m n_i-m\]

Tukey-Kramer: distintos tamaños muestrales

\[q_N = \frac{X_{i.} - X_{j.}}{\sqrt{MSE \left(\frac{\frac{1}{n_i}+ \frac{1}{n_j}}{2}\right)}}, \qquad df =N = \sum_{i=1}^m n_i - m\]

En ambos casos \(q_N\) sigue una distribución de rangos student (studendized range distribution) que es similar a t-student (studentized t-test) pero diferente. Usabmos el \(\alpha\) original (e.g., 0.05). Segimos el proceso de de valor crítico o p-value.

El Tukey es menos conservador que Bonferroni. Tiene un control estricto sobre error de Tipo I (menos estricto que Bonferroni). Tiene buena potencia.

14.5. Análisis de asociación en v.a. discretas

En esta sección se presentan algunos tests que permiten analizar la asociación entre variables aleatorias discretas.

14.5.1. Tabla de contingencia

Consideremos que tenemos observaciones de dos variables discretas o categóricas (o bien variables continuas agrupadas en clases) \(X\) y \(Y\) referidas a una muestra. La tabla de contingencia contiene los efectivos correspondientes a cada pareja de valores (o clases) de ambas variables:

\[\begin{split}\begin{array}{|c||c|c|c|c|c||c|} \hline \bf{X\mid Y}& \bf{d_1} &\cdots &\bf{d_k}& \cdots&\bf{d_s}& \text{total}\\ \hline\hline \bf{c_1} & n_{11} & \cdots & n_{1k}& \cdots& n_{1s} &\bf { n_{1\bullet}}\\ \hline \vdots & \vdots &&\vdots&&\vdots&\vdots\\ \hline \bf{c_h} & n_{h1} & \cdots & n_{hk}& \cdots& n_{hs} & \bf{n_{h\bullet}}\\ \hline \vdots & \vdots &&\vdots&&\vdots&\vdots\\ \hline \bf{c_r} & n_{r1} & \cdots & n_{rk}& \cdots& n_{rs} & \bf{n_{r\bullet}}\\ \hline \hline \text{total} &\bf{n_{\bullet 1}}&\cdots& \bf{n_{\bullet k}}&\cdots&\bf{n_{\bullet s}}&\bf{n}\\ \hline \end{array}\end{split}\]

Cada linea y cada columna corresponden a una submuestra particular. La fila de índice \(h\) es la distribución sobre \( d_1,\ldots,d_s\) de los individuos para los cuales la variable \(X\) toma el valor \(c_h\). La columna de índice \(k\) es la distribución sobre \(c_1,\ldots,c_r\) de los individuos para los cuales las variable \( Y\) toma el valor \(d_k\).

Los totales de renglones y columnas en la tabla \(n_{h\bullet}, n_{\bullet k}\) (etc.) se denominan frecuencias marginales.

Los valores en las celdas internas \(n_{hk}\) (etc.) se denominan frecuencias observadas.

14.5.2. Test de Chi-cuadrado de independencia para tablas de contingencia

Tenemos observaciones de dos variables discretas o categóricas (o bien variables continuas agrupadas en clases) \(X\) y \(Y\) referidas a una muestra de distribuciones desconocidas. Las observaciones de \(X\) son independendes, y las observaciones de \(Y\) son independendes. Nos interesa si \(X\) y \(Y\) son independientes. Tenemos las hipótesis:

\[\begin{split}H_0: \forall h \in [1,r], k \in [1, s], \qquad P_{hk} = p_h q_k \\ H_1: \exists h \in [1,r], k \in [1, s], \qquad P_{hk} \neq p_h q_k\end{split}\]

donde

• \(P_{hk}\) se refiere a la probabilidad conjunta de que \(X\) sea \(c_h\) y \(Y\) sea \(d_k\)

• \(p_{h}\) se refiere a la probabilidad marginal de que \(X\) sea \(c_h\)

• \(q_{k}\) se refiere a la probabilidad marginal de que \(Y\) sea \(d_k\)

Ya que no conocemos sus distribuciones, estimamos \(P_{hk}, p_{h}, q_{k}\) a partir de la muestra, utilizando la tabla de contingencia:

\[ \hat P_{hk} = \frac{n_{hk}}{n} \qquad \hat p_{h} = \frac{n_{h \bullet}}{n} \qquad \hat q_{k} = \frac{n_{\bullet k}}{n} \]

(Por qué podemos estimar así?)

Asi que bajo \(H_0\), queremos que los valores de \(\hat P_{hk}\) y \(\hat p_{h} \hat q_{k}\) estén cercas:

\[ \frac{n_{hk}}{n} \approx \frac{n_{h \bullet}}{n} \frac{n_{\bullet k}}{n} \qquad \Leftrightarrow \qquad n_{hk} \approx \frac{n_{h \bullet} n_{\bullet k}}{n} \]

Recuerda que \(n_{hk}\) es la frecuencia observada (\(O_{hk}\)), y aquí el producto \(\displaystyle \frac{n_{h \bullet} n_{\bullet k}}{n}\) se denomina la frecuencia esperada (\(E_{hk}\)). Definimos un estadístico (v.a.) que representa la distancia entre ellos para todas combinaciones de valores:

\[TS = n \sum\limits_{k=1}^s \sum_{h=1}^r \frac{(\hat P_{hk} - \hat p_{h} \hat q_{k})^2}{\hat p_{h} \hat q_{k}} = \sum\limits_{k=1}^s \sum_{h=1}^r \frac{(O_{hk} - E_{hk})^2}{E_{hk}} = \sum\limits_{k=1}^s \sum_{h=1}^r \frac{\left(n_{hk} - \displaystyle \frac{n_{h \bullet} n_{\bullet k}}{n}\right)^2}{\displaystyle \frac{n_{h \bullet} n_{\bullet k}}{n}} \]

Tenemos un teorema que dice bajo \(H_0\), este estadístico \(TS\) converge hacia una distribución \(\chi^2_{(r-1)(s-1)}\) cuando \( n\) tiende al infinito.

Por qué el grado de libertad es \((r-1)(s-1)\)?

Ejemplo (campus y tiempo libre)

Consideremos dos variables categóricas referidas a una muestra de los estudiantes de la UACh, una de las cuales indica el campus, y la otra indica que hace en su tiempo libre. Se trata de analizar si existe una dependencia entre el campus y actividades en su tiempo libre. Supongamos que la tabla de contingencia observada es la siguiente (los números son inventados):

\[\begin{split}\begin{array}{|c||c|c|c|c|} \hline \bf{X\mid Y}& \text{Deporte} & \text{Estudio} &\text{Trabajo} & \text{Total} \\ \hline\hline \text{Campus Isla Teja} & 68 & 56 & 32 & 156\\ \hline \text{Campus Miraflores} & 52 & 72 & 20 & 144\\ \hline \text{Total} &120 & 128 & 52 & 300\\ \hline \end{array}\end{split}\]

En tu tarea, tienes que desarrollar el test. Aquí, mostramos el resultado de R:

#usando el test predifinido en R
table <- rbind(c(68, 56, 32), c(52, 72, 20))
chisq.test(table, correct=FALSE)

	Pearson's Chi-squared test

data:  table
X-squared = 6.4, df = 2, p-value = 0.04

¿Cómo calculamos el p-value arriba con la función pchisq en R?

¿Cómo calculamos el valor crítico (\(\alpha=0.05\)) con la función qchisq en R?

14.5.3. Test de Bowker (o McNemar-Bowker)

Este test se utiliza para analizar la simetría en una tabla de contingencia construida a partir de dos muestras dependientes (\(X, Y\)), típicamente pre y post un tratamiento.

Las hípotesis de este test es:

\[\begin{split}H_0: \forall i, j \in [1,k], i \neq j, \qquad p_{ij} = p_{ji} \\ H_1: \exists i, j \in [1,k], i \neq j, \qquad p_{ij} \neq p_{ji}\end{split}\]

donde

• \(p_{ij}\) se refiere a la probabilidad de que \(X\) sea \(i\) y \(Y\) sea \(j\) (usamos \(i\) o \(j\) para referir a un valor de una clase o categoría, aunque strictamente debemos escribir \(c_i\) o \(c_j\))

• \(p_{ji}\) se refiere a la probabilidad de que \(X\) sea \(j\) y \(Y\) sea \(i\)

El estadístico del test de Bowker se construye como:

\[ B= \sum_{j>i} \frac{(n_{ij} - n_{ji})^2}{n_{ij} + n_{ji}}\]

De manera que bajo la hipotesis nula de simetría en la tabla de contingencia, es decir que no hay efecto del tratamiento, se cumple asintóticamente que:

\[B \sim \chi^2_{\frac{k(k-1)}{2}}\]

Si se rechaza la hipótesis nula (cuándo?), y además se observan mas efectivos sobre la diagonal (upper off-diagonal) que bajo la diagonal (lower off-diagonal) se puede concluir que hay un efecto positivo del tratamiento en el criterio en análisis.

Ejemplo:

Considere los datos de 54 estudiantes sometidos a una evaluación sobre sus habilidad de cohesión en la redacción de textos, antes y después de una intervención didáctica. La evaluación utiliza una escala de 4 niveles desde no logrado a totalmente logrado. Se pide evaluar si la intervención ha tenido efecto o no en la habilidad estudiada.

tabla <- rbind(c(0,2,3,2), c(2,2,3,7), c(3,4,6,19), c(0,0,1,2))
print(tabla)
B = 0
efsd = efbd = 0
for (i in 1:3){
    for (j in (i+1):4){
        B = B + (tabla[i,j] - tabla[j,i])^2 / (tabla[i,j] + tabla[j,i])
        efsd = efsd + tabla[i,j]
        efbd = efbd + tabla[j,i]
    }
}
cat("B:", B, "\n")
cat("sobre diag vs. bajo diag:", c(efsd, efbd), "\n")

k = 4*3/2
q = qchisq(0.95, k)
print(q)

p <- 1-pchisq(B, k)
print(p)

     [,1] [,2] [,3] [,4]
[1,]    0    2    3    2
[2,]    2    2    3    7
[3,]    3    4    6   19
[4,]    0    0    1    2
B: 25.34286 
sobre diag vs. bajo diag: 36 10 
[1] 12.59159
[1] 0.000294974

Se rechaza la hipótesis nula de simetría en favor de un efecto positivo de la intervención didáctica.

##usando el test predefinido en R
mcnemar.test(tabla)

	McNemar's Chi-squared test

data:  tabla
McNemar's chi-squared = 25.343, df = 6, p-value = 0.000295